Simulating Tetrathionate Circuit with Bioscrape and BioCRNpyler¶
Liana Merk¶
BE240 Spring 2020
In [1]:
from bioscrape.types import Model
from bioscrape.sbmlutil import import_sbml
from bioscrape.simulator import py_simulate_model
from biocrnpyler import *
import numpy as np
import pylab as plt
import tqdm
import pandas as pd
import matplotlib as mpl
import holoviews as hv
import bokeh.plotting
import bokeh.io
from bokeh.models import Range1d
from bokeh.layouts import gridplot
bokeh.io.output_notebook()
hv.extension('bokeh')
import seaborn as sns
sns.set()
BioCRNpyler Custom Classes¶
Phosphorylation¶
In [2]:
class Phosphorylation(Mechanism):
def __init__(self,name="phosphorylation",mechanism_type="phosphorylation"):
Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
def update_species(self,s1, **keywords):
phosphoprotein = ComplexSpecies([s1,Species("P",material_type="phosphate")])
return [phosphoprotein,Species("P",material_type="phosphate")]
class AutoPhosphorylation(Phosphorylation):
def __init__(self,name="autophosphorylation",mechanism_type="phosphorylation"):
Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
def update_reactions(self,s1,component = None, kphos = None, kdephos = None, \
part_id = None,**keywords):
if part_id == None:
part_id = s1.name+"-P"
if kphos == None and component != None:
try:
kphos = component.get_parameter("kphos",part_id = part_id,mechanism = self)
except ValueError:
kphos = 0
if kdephos == None and component != None:
try:
kdephos = component.get_parameter("kdephos",part_id = part_id,mechanism = self)
except ValueError:
kdephos = 0
#if component == None and (kphos == None or kdephos == None):
# raise ValueError("Must pass in a Component or values for kphos, kdephos.")
if((kphos == 0 and kdephos == 0) or (kphos == None and kdephos == None)):
print(component)
print("kphos is {} and kdephos is {}".format(kphos,kdephos))
raise ValueError("Must pass in at least one value, kphos or kdephos for "+part_id)
rxns = []
phosphate = Species("P",material_type="phosphate")
phosphoprotein = self.update_species(s1)[0]
if(kphos == None):
kphos = 0
if(kdephos == None):
kdephos = 0
if(kdephos > 0 and kphos > 0):
rxns += [Reaction([s1,phosphate],[phosphoprotein],k = kphos,k_rev=kdephos)]
elif(kphos > 0):
rxns += [Reaction([s1,phosphate],[phosphoprotein],k = kphos)]
elif(kdephos > 0):
rxns += [Reaction([phosphoprotein],[s1,phosphate],k = kdephos)]
return rxns
class SimpleKinaseMechanism(Phosphorylation):
def __init__(self,name="Kinase",mechanism_type="phosphorylation"):
Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
self.phosphate = Species("P",material_type="phosphate")
def update_species(self,s1,s2):
#s1 causes s2 to get phosphorylated. So, we need to update species
#with s2's phosphorylated form
return [ComplexSpecies([s2,self.phosphate]),self.phosphate]
def update_reactions(self,s1,s2,component = None, k = None, \
part_id = None,**keywords):
if part_id == None:
part_id = s1.name+"-"+s2.name
if k == None and component != None:
k = component.get_parameter("kkinase",part_id = part_id,mechanism = self)
s2_p = ComplexSpecies([s2,self.phosphate])
return [Reaction([s1,s2],[s1,s2_p],k=k)]
class TwoStepKinaseMechanism(Phosphorylation):
def __init__(self,name="Kinase",mechanism_type="phosphorylation"):
Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
self.phosphate = Species("P",material_type="phosphate")
def update_species(self,s1,s2):
#s1 causes s2 to get phosphorylated. So, we need to update species
#with s2's phosphorylated form
s2_p = ComplexSpecies(s2,self.phosphate)
s1_s2_complex = ComplexSpecies([s1,s2])
s1_s2_plus_p = ComplexSpecies([s1,s2_p])
return [s2_p,s1_s2_complex,s1_s2_plus_p,self.phosphate]
def update_reactions(self,s1,s2,component = None, kb1 = None, ku1=None,kkinase=None, kb2 = None, ku2=None, \
part_id = None,**keywords):
"""this mechanism involves three reactions:
1) s1 binds to s2 to make Complex(s1,s2) this uses kb1 and ku1 for binding and unbinding
2) Complex(s1,s2) + phosphate makes Complex(s1,Complex(s2,phosphate)) this uses kkinase, and only goes forwards
3) Complex(s1,Complex(s2,phosphate)) unbinds to form s1, Complex(s2,phosphate). this uses kb2 and ku2
"""
if part_id == None:
part_id = s1.name+"-P-"+s2.name
if kkinase == None and component != None:
kkinase = component.get_parameter("kkinase",part_id = part_id,mechanism = self)
if ku1 == None and component != None:
ku1 = component.get_parameter("ku1",part_id = part_id,mechanism = self)
if kb1 == None and component != None:
kb1 = component.get_parameter("kb1",part_id = part_id,mechanism = self)
if kb2 == None and component != None:
kb2 = component.get_parameter("kb2",part_id = part_id,mechanism = self)
if ku2 == None and component != None:
ku2 = component.get_parameter("ku2",part_id = part_id,mechanism = self)
specs = self.update_species(s1,s2)
binding = Reaction([s1,s2],specs[1],k=kb1,k_rev=ku1)
kinase = Reaction([specs[1],self.phosphate],[specs[2]],k=kkinase)
kunbinding = Reaction([specs[2]],[s1,specs[0]],k=ku2,k_rev=kb2)
return [binding,kinase,kunbinding]
class PhosphoTransfer(Phosphorylation):
def __init__(self,name="phosphotransfer",mechanism_type="phosphorylation"):
Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
def update_reactions(self,s1,s2,component = None, k = None, \
part_id = None,**keywords):
if part_id == None:
part_id = s1.name+"-"+s2.name
if k == None and component != None:
k = component.get_parameter("ktransfer",part_id = part_id,mechanism = self)
s1_p = self.update_species(s1)[0]
s2_p = self.update_species(s2)[0]
return [Reaction([s1_p,s2],[s2_p,s1],k=k)]
Combinatorial Binding¶
In [3]:
class Combinatorial_Cooperative_Binding(Mechanism):
"""a reaction where some number of binders bind combinatorially to a bindee"""
def __init__(self,name="Combinatorial_Cooperative_binding",
mechanism_type="cooperative_binding"):
Mechanism.__init__(self,name,mechanism_type)
def make_cooperative_complex(self,combo,bindee,cooperativity):
"""given a list of binders and their cooperativities, make a complex
that contains the binders present N number of times where N is
each one's cooperativity"""
complexed_species_list = []
for binder in combo:
binder_cooperativity = int(cooperativity[binder.name])
#I hope that cooperativity is an int! what if it isn't
if(binder_cooperativity > 1):
complexed_species_list += [Multimer(binder,binder_cooperativity)]
else:
complexed_species_list += [binder]
complexed_species_list += [bindee]
if(len(complexed_species_list)==1):
myspecies = complexed_species_list[0]
else:
myspecies = ComplexSpecies(complexed_species_list)
return myspecies
def update_species(self,binders,bindee,cooperativity=None,\
component = None, part_id = None, **kwords):
cooperativity_dict = {}
for binder in binders:
binder_partid = part_id+"_"+binder.name
if ((cooperativity == None) or (type(cooperativity)==dict and binder_partid not in cooperativity) \
and (component != None)):
coop_val = component.get_parameter("cooperativity", part_id = binder_partid, mechanism = self)
elif type(cooperativity)==dict and binder_partid in cooperativity:
coop_val = cooperativity[binder_partid]
if component == None and ( cooperativity == None):
raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
cooperativity_dict[binder.name]=coop_val
#allbinders = binders+[bindee]
out_species = []
for i in range(1, len(binders)+1):
for combo in it.combinations(binders,i):
#go through every possible combination of reactants and dna and make
#all the complexes
out_species += [self.make_cooperative_complex(combo,bindee,cooperativity_dict)]
return out_species
def update_reactions(self,binders,bindee,component=None,kbs=None,kus=None,\
part_id = None,cooperativity=None,**kwords):
binder_params = {}
for binder in binders:
binder_partid = part_id+"_"+binder.name
if ((type(kbs)==dict and binder not in kbs) or (type(kbs)!=dict and component != None)):
kb = component.get_parameter("kb", part_id = binder_partid, mechanism = self)
elif(type(kbs)==dict and binder in kbs):
kb = kbs[binder.name]
elif(type(kbs)!=dict and component == None):
raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
if ((type(kus)==dict and binder not in kus) or (kus == None and component != None)):
ku = component.get_parameter("ku", part_id = binder_partid, mechanism = self)
elif(type(kus)==dict and binder in kus):
ku = kus[binder.name]
elif(type(kus)!=dict and component == None):
raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
if ((cooperativity == None) or (type(cooperativity)==dict and binder.name not in cooperativity) \
and component != None):
coop_val = component.get_parameter("cooperativity", part_id = binder_partid, mechanism = self)
elif type(cooperativity)==dict and binder.name in cooperativity:
coop_val = cooperativity[binder.name]
if component == None and (kb == None or ku == None or cooperativity == None):
raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
binder_params[binder] = {"kb":kb,"ku":ku,"cooperativity":coop_val}
#out_rxns = []
rxndict = {}
coop_dict = {a.name:binder_params[a]["cooperativity"] for a in binder_params}
for i in range(1, len(binders)+1):
for combo in it.combinations(binders,i):
#come up with all combinations of binders
product = self.make_cooperative_complex(combo,bindee,coop_dict)
#this is the complex which becomes the product
for binder in combo:
#then in each case, one binder can be added to make
#this combination
reactant = tuple(set(combo)-set([binder]))
rxn_prototype = (binder,reactant)
#this part makes a describer of the reaction; which reactants are combining?
#print(rxn_prototype)
#print(rxndict)
if(rxn_prototype in rxndict):
#if we already did this reaction then forget about it
continue
else:
reactant_complex = self.make_cooperative_complex(reactant,bindee,coop_dict)
reaction = Reaction(inputs=[binder, reactant_complex], outputs=[product],
input_coefs=[binder_params[binder]["cooperativity"], 1], output_coefs=[1], \
k=binder_params[binder]["kb"],k_rev=binder_params[binder]["ku"])
rxndict[rxn_prototype]=reaction
return [rxndict[a] for a in rxndict]
In [4]:
class CombinatorialPromoter(Promoter):
def __init__(self, name, regulators, leak = True, assembly = None,
transcript = None, length = 0, mechanisms = {},
parameters = {},tx_capable_list = None,cooperativity = None, **keywords):
"""
tx_capable_list = [[1,2,3],[0,2,3]] this means having regulator 1, 2, and 3 will transcribe
but also 3, 2, and 0.
#TODO make it force sorted list
"""
if not isinstance(regulators, list):
#you could give one string as a regulator
regulators = [regulators]
self.cooperativity = cooperativity
self.regulators = []
for regulator in regulators:
if(isinstance(regulator,str)):
self.regulators += [self.set_species(regulator, material_type = "protein")]
#if it's a string then assume it's a protein
elif(isinstance(regulator,Species)):
#if it's already a species then add it wholesale
self.regulators += [regulator]
#after we've sanitized the inputs, then sort
self.regulators = sorted(self.regulators)
#now let's work out the tx_capable_list
if(tx_capable_list == None):
#if nothing is passed assume default
self.tx_capable_list = [set([a.name for a in self.regulators])]
elif(type(tx_capable_list)==list):
#if the user passed a list then the user knows what they want
self.tx_capable_list = [set(a) for a in tx_capable_list]
self.leak = leak
self.default_mechanisms = {"binding": Combinatorial_Cooperative_Binding()}
Promoter.__init__(self, name = name, assembly = assembly,
transcript = transcript, length = length,
mechanisms = mechanisms, parameters = parameters,
**keywords)
self.complex_combinations = {}
self.tx_capable_complexes = []
def update_species(self):
mech_tx = self.mechanisms["transcription"]
mech_b = self.mechanisms['binding']
#set the tx_capable_complexes to nothing because we havent updated species yet!
self.tx_capable_complexes = []
species = []
self.complexes = []
if self.leak is not False:
species += mech_tx.update_species(dna = self.assembly.dna)
bound_species = mech_b.update_species(self.regulators,self.assembly.dna,\
component = self,part_id = self.name,cooperativity=self.cooperativity)
#above is all the species with DNA bound to regulators. Now, we need to extract only the ones which
#are transcribable
for bound_complex in bound_species:
species_inside = []
for regulator in self.regulators:
if(regulator.name in bound_complex.name):
species_inside += [regulator.name]
#print("species_inside is "+str(species_inside))
#print("tx capable is " +str(self.tx_capable_list))
#print(set(species_inside) in [set(a) for a in self.tx_capable_list])
if(set(species_inside) in [set(a) for a in self.tx_capable_list]):
tx_capable_species = mech_tx.update_species(dna = bound_complex, transcript = self.transcript, \
protein = self.assembly.protein)
species +=tx_capable_species[1:]
self.tx_capable_complexes +=[bound_complex]
species+=bound_species
return species
def update_reactions(self):
reactions = []
mech_tx = self.mechanisms["transcription"]
mech_b = self.mechanisms['binding']
if self.leak != False:
reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
reactions += mech_b.update_reactions(self.regulators,self.assembly.dna,component = self,\
part_id = self.name,cooperativity=self.cooperativity)
if(self.tx_capable_complexes == None or self.tx_capable_complexes == []):
#this could mean we haven't run update_species() yet
species = self.update_species()
if(self.tx_capable_complexes == []):
#if it's still zero after running update_species then we could be in trouble
warn("nothing can transcribe from combinatorial promoter {}".format(self.name))
else:
for specie in self.tx_capable_complexes:
tx_partid = self.name
for part in specie.species_set:
#construct the name of the promoter with regulators bound
if part.material_type == "dna":
#the DNA doesn't matter
pass
else:
#put in the regulators!
tx_partid += "_"+part.name
if(tx_partid[0]=="_"):
tx_partid = tx_partid[1:]
#if it's bound to RNAP then it transcribes, right?
tx_partid = tx_partid+"_RNAP"
#print("tx_partid is "+ str(tx_partid))
reactions += mech_tx.update_reactions(dna = specie, component = self, part_id = tx_partid, \
transcript = self.transcript, protein = self.assembly.protein)
return reactions
Tetrathionate Sensing Circuit in Bioscrape¶
Sensory kinase, ttrS, is membrane bound and will phosphorylate when tetrathionate is present. In turn, it phosphorylates ttrR. A ttrR-phosphorylated dimer is needed to activate the pTTR promoter.
Let's create a mass action model for the system.
\begin{align}
\textrm{Polymerase binds to P1d: } P + P1d &\underset{k^u_1}{\overset{k^b_1}\rightleftharpoons} P1d : P \\
\textrm{Transcription of ttrS, ttrR: } P1d : P &\xrightarrow{k_{tx}} ttrS_T + ttrR_T + P1d + P\\
\textrm{Translation ttrS_T to ttrS: } ttrS_T + R &\underset{k^u_2}{\overset{k^b_2}\rightleftharpoons} ttrS_T: R \xrightarrow{k_{tl}} ttrS_T + R + ttrS \\
\textrm{Translation ttrR_T to ttrR: } ttrR_T + R &\underset{k^u_3}{\overset{k^b_3}\rightleftharpoons} ttrR_T: R \xrightarrow{k_{tl}} ttrR_T + R + ttrR \\
\textrm{ttrS phosphorylation: } ttrS + tt &\underset{k^u_4}{\overset{k^b_4}\rightleftharpoons} ttrS : P + tt\\
\textrm{ttrR binding to unphosphorylated ttrS: } ttrR + ttrS &\underset{k^u_5}{\overset{k^b_5}\rightleftharpoons} ttrR:ttrS\\
\textrm{ttrR binding to phosphorylated ttrS: } ttrR + ttrS:P &\underset{k^u_6}{\overset{k^b_6}\rightleftharpoons} ttrR:ttrS:P\\
\textrm{Phosphorylation of ttrR: } ttrR:ttrS:P &\underset{k^u_7}{\overset{k^b_7}\rightleftharpoons} ttrR:P + ttrS\\
\textrm{Dephosphorylation of ttrR: } ttrR:P &\xrightarrow{k_{dephos}} ttrR \\
\textrm{Dimerization of ttrR:P: } 2ttrR:P &\underset{k^u_8}{\overset{k^b_8}\rightleftharpoons} ttrR:P_{dimer} \\
\textrm{ttrR:P dimer binding to pttr promoter: } ttrR:P_{dimer} + pttr0 &\underset{k^u_9}{\overset{k^b_9}\rightleftharpoons} pttr1 \\
\textrm{Polymerase binds to pttr: } P + pttr1 &\underset{k^u_{10}}{\overset{k^b_{10}}\rightleftharpoons} pttr1 : P \\
\textrm{Transcription of GFP: } pttr1 : P &\xrightarrow{k_{tx}} pttr1 + P + GFP_T \\
\textrm{Translation GFP: } GFP_T + R &\underset{k^u_{11}}{\overset{k^b_{11}}\rightleftharpoons} GFP_T: R \xrightarrow{k_{tl}} GFP_T + R + GFP \\
ttrS_T &\xrightarrow{\delta} \emptyset\\
ttrS &\xrightarrow{\delta} \emptyset\\
ttrR_T &\xrightarrow{\delta} \emptyset\\
ttrR &\xrightarrow{\delta} \emptyset\\
GFP_T &\xrightarrow{\delta} \emptyset \\
GFP &\xrightarrow{\delta} \emptyset
\end{align}
In [5]:
#Species: A list of all the species
species = ["P", "P1d", "P1d:P", "ttrS_T", "ttrR_T", "R", "ttrR_T:R", "ttrS_T:R", "tt", "ttrR", "ttrS",
"ttrS:P", "ttrR:ttrS", "ttrR:ttrS:P", "ttrR:P", "ttrR:P_dimer", "pttr0", "pttr1", "pttr1:P", "GFP_T", "GFP_T:R", "GFP"]
In [6]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 50), ("k2u", 10),
("k3b", 30), ("k3u", 10),
("k4b", 80), ("k4u", 2),
("k5b", 50), ("k5u",1),
("k6b", 50), ("k6u", 1),
("k7b", 50), ("k7u", 1),
("k8b", 80), ("k8u", 1),
("k9b", 50), ("k9u", 0),
("k10b", 50), ("k10u",1),
("k11b", 10), ("k11u", 1),
("ktx", .5), ("ktl", 5), ("k_dephos", 0.1), ("delta", .5)]
In [7]:
#Reactions: A list of reactions [rxn1, rxn2...]. Each reaction is a tuple ([Input Species], [Output Species], "propensity_type", {propensity_parameters})
# P + P1d <-> P1d:P
rxn1b, rxn1u = (["P", "P1d"], ["P1d:P"], "massaction", {"k":"k1b"}), (["P1d:P"], ["P", "P1d"], "massaction", {"k":"k1u"})
# P1d:P -> ttrS_T + ttrR_T
rxntx1 = (["P1d:P"], ["ttrS_T", "ttrR_T"], "massaction", {"k":"ktx"})
# ttrS_T + R <-> ttrS_T:R
rxn2b, rxn2u = (["ttrS_T", "R"], ["ttrS_T:R"], "massaction", {"k":"k2b"}), (["ttrS_T:R"], ["ttrS_T", "R"], "massaction", {"k":"k2u"})
# ttrS_T:R -> ttrS_T + R + ttrS
rxntl1 = (["ttrS_T:R"], ["ttrS_T", "ttrS", "R"], "massaction", {"k":"ktl"})
# ttrR_T + R <-> ttrR_T:R
rxn3b, rxn3u = (["ttrR_T", "R"], ["ttrR_T:R"], "massaction", {"k":"k3b"}), (["ttrR_T:R"], ["ttrR_T", "R"], "massaction", {"k":"k3u"})
# ttrR_T:R -> ttrR_T + R + ttrR
rxntl2 = (["ttrR_T:R"], ["ttrR_T", "ttrR", "R"], "massaction", {"k":"ktl"})
# ttrS + tt <-> ttrS:P + tt
rxn4b, rxn4u = (["ttrS", "tt"], ["ttrS:P", "tt"], "massaction", {"k":"k4b"}), (["ttrS:P", "tt"], ["ttrS", "tt"], "massaction", {"k":"k4u"})
# ttrR + ttrS <-> ttrR:ttrS
rxn5b, rxn5u = (["ttrR", "ttrS"], ["ttrR:ttrS"], "massaction", {"k":"k5b"}), (["ttrR:ttrS"], ["ttrR", "ttrS"], "massaction", {"k":"k5u"})
# ttrR + ttrS <-> ttrR:ttrS:P
rxn6b, rxn6u = (["ttrS:P", "ttrR"], ["ttrR:ttrS:P"], "massaction", {"k":"k6b"}), (["ttrR:ttrS:P"], ["ttrS:P", "ttrR"], "massaction", {"k":"k6u"})
# ttrR:ttrS:P <-> ttrR:P + ttrS
rxn7b, rxn7u = (["ttrR:ttrS:P"], ["ttrR:P", "ttrS"], "massaction", {"k":"k7b"}), (["ttrR:P", "ttrS"], ["ttrR:ttrS:P"], "massaction", {"k":"k7u"})
# ttrR:P -> ttrR
rxn_dephos = (["ttrR:P"], ["ttrR"], "massaction", {"k":"k_dephos"})
# 2ttrR:P <-> ttrR:P_dimer
rxn8b, rxn8u = (["ttrR:P", "ttrR:P"], ["ttrR:P_dimer"], "massaction", {"k":"k8b"}), (["ttrR:P_dimer"], ["ttrR:P", "ttrR:P"], "massaction", {"k":"k8u"})
# ttrR:P_dimer + pttr0 <-> pttr1
rxn9b, rxn9u = (["ttrR:P_dimer", "pttr0"], ["pttr1"], "massaction", {"k":"k9b"}), (["pttr1"], ["ttrR:P_dimer", "pttr0"], "massaction", {"k":"k9u"})
# P + pttr1 <-> pttr1:P
rxn10b, rxn10u = (["P", "pttr1"], ["pttr1:P"], "massaction", {"k":"k10b"}), (["pttr1:P"], ["P", "pttr1"], "massaction", {"k":"k10u"})
# pttr1:P -> pttr1 + P + GFP_T
rxntx2 = (["pttr1:P"], ["pttr1", "P", "GFP_T"], "massaction", {"k":"ktx"})
# GFP_T + R <-> GFP_T:R
rxn11b, rxn11u = (["GFP_T", "R"], ["GFP_T:R"], "massaction", {"k":"k11b"}), (["GFP_T:R"], ["GFP_T", "R"], "massaction", {"k":"k11u"})
# GFP_T:R -> GFP_T + GFP + R
rxntl3 = (["GFP_T:R"], ["GFP_T", "GFP", "R"], "massaction", {"k":"ktl"})
# RNA deg
rxnd_ttrS_T = (["ttrS_T"], [], "massaction", {"k":"delta"})
rxnd_ttrR_T = (["ttrR_T"], [], "massaction", {"k":"delta"})
rxnd_GFP_T = (["GFP_T"], [], "massaction", {"k":"delta"})
# Protein Deg
rxnd_ttrS = (["ttrS"], [], "massaction", {"k":"delta"})
rxnd_ttrR = (["ttrR"], [], "massaction", {"k":"delta"})
rxnd_GFP = (["GFP"], [], "massaction", {"k":"delta"})
reactions = [rxn1b, rxn1u, rxn2b, rxn2u, rxn3b, rxn3u, rxn4b, rxn4u,
rxn5b, rxn5u, rxn6b, rxn6u, rxn7b, rxn7u, rxn8b, rxn8u,
rxn9b, rxn9u, rxn10b, rxn10u, rxn11b, rxn11u,
rxntx1, rxntx2,
rxntl1, rxntl2, rxntl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
In [8]:
#An initial condition for each species (uninitialized species default to 0)
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":25,
"R":100,
}
In [9]:
#Instantiate the Model [The only object must of us will care about]
M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
timepoints = np.linspace(0, 600, 100)
Results = py_simulate_model(timepoints, Model = M, stochastic = False) #py_simulate_model takes care of all other objects for you
Results.head()
Out[9]:
Matplotlib plotting of Bioscrape results¶
In [10]:
#Bioscrape returns a Pandas Dataframe by default, so plotting is easy!
plt.figure(figsize = (20, 10))
plt.subplot(221)
plt.title("Tetrathionate Regulators")
plt.plot(timepoints, Results["ttrR"], color = "pink", label = "ttrR")
plt.plot(timepoints, Results["ttrR_T"] + Results["ttrR_T:R"],":", color = "pink", label = "ttrR_T")
plt.plot(timepoints, Results["ttrS"], color = "blue", label = "ttrS")
plt.plot(timepoints, Results["ttrS_T"] + Results["ttrS_T:R"], ":", color = "blue",label = "ttrS_T")
plt.plot(timepoints, Results["tt"], color = "red", label = "tt")
plt.legend()
plt.subplot(222)
plt.title("Phosphorylation Pathway")
plt.plot(timepoints, Results["ttrR:ttrS"], label = "ttrR:ttrS")
plt.plot(timepoints, Results["ttrR:ttrS:P"], label = "ttrR:ttrS:P")
plt.plot(timepoints, Results["ttrR:P"], label = "ttrR:P")
plt.plot(timepoints, Results["ttrR:P_dimer"], label = "ttrR:P dimer")
plt.legend()
plt.subplot(223)
plt.title("Promoters")
plt.plot(timepoints, Results["P1d:P"], color = 'blue', label = "P1d:P")
plt.plot(timepoints, Results["P1d"], ":", color = 'blue', label = "P1d")
plt.plot(timepoints, Results["pttr1"], color = "cyan", label = "pttr1")
plt.plot(timepoints, Results["pttr0"], ":", color = "cyan", label = "pttr0")
plt.legend()
plt.subplot(224)
plt.title("Readout")
plt.plot(timepoints, Results["GFP_T"] + Results["GFP_T:R"], ":", color = 'green', label = "GFP mRNA")
plt.plot(timepoints, Results["GFP"], color = 'green', label = "GFP")
plt.legend()
Out[10]:
Holoviews plotting of Bioscrape results¶
Sort out Plotting¶
In [11]:
def plot_holoviews_ttr_bioscrape(results_melted):
plot_group = []
for row in results_melted.iterrows():
sp = row[1]['species']
if sp in ["P1d", "P1d:P", "pttr0", "pttr1", "pttr1:P"]:
plot_group.append('promoters')
elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T:R"]:
plot_group.append('mRNA')
elif sp in ["ttrS:P", "ttrR:ttrS", "ttrR:ttrS:P", "ttrR:P", "ttrR:P_dimer"]:
plot_group.append('complexes')
elif sp in ["ttrR", "ttrS", "GFP"]:
plot_group.append('proteins')
elif sp in ["P", "R", "tt"]:
plot_group.append('machinery and signals')
else:
print(sp)
results_melted['plot_group'] = plot_group
plot_list = []
for group in results_melted['plot_group'].unique():
plot_list.append(hv.Points(
data=results_melted.loc[results_melted['plot_group'] == group],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
[ 'species']
).opts(
tools=['hover'],
legend_position = 'right',
title = f'{group}',
width = 500,
#shared_axes=False
).overlay('species'))
return plot_list
Now, proceed with plotting¶
In [12]:
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
Out[12]:
In [13]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(5).opts(shared_axes=False)
Out[13]:
Decrease ttrS phosphorylation¶
Now, I am curious how the time course of the complexes changes if we decrease the sensory kinase's rate of phosphorylating. Perhaps the cell is low on phosphates due to other processes, and thus this rate is lower. We will thus just change the parameter associated with the forward reaction of ttrS phosphorylation.
In [14]:
M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
# See what happens if our sensory molecule is slower
M.set_parameter("k4b",0.0001)
#Simulate the Model
timepoints = np.linspace(0, 600, 100)
Results = py_simulate_model(timepoints, Model = M, stochastic = False)
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
We can use our plotting function we defined above:
In [15]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(5).opts(shared_axes=False)
Out[15]:
As expected, the entire response involves a delay now. Since phosphorylation is slower, we actually have no accumulation of ttrS:P. It is used immediately after it is created, and never builds up. We do have a lot of ttrR:ttrS compelxes, which are "stuck"
Test if Dose-Response Curve agrees with experimental data¶
In [16]:
M.set_parameter("delta", .5)
GFP_max = []
#Different initial values of S
tts = np.logspace(-3, 3, 100)
for tt0 in tts:
x0["tt"] = tt0 #Change my initial condition dictionary
M.set_species(x0)
Results = py_simulate_model(timepoints, Model = M)
GFP_max.append(np.max(Results["GFP"]))
In [17]:
LM_19 = pd.read_csv('LM19_data.csv')
LM_19
Out[17]:
Due to the unknown calibration factor between fluorescence and concentration, I have applied transformations to get the data in what seems to qualitatively match this curve. This is preliminary, and I hope to find a better way to show this agreement between data and simulation.
In [18]:
p = bokeh.plotting.figure(
height=300,
width=550,
x_axis_label='[tt]',
y_axis_label='Max GFP',
x_axis_type = 'log',
title = 'Comparison of Experimental Sigmoid with Simulated'
)
p.line(x=tts,y=GFP_max,)
# Try to transform data in some way that we can compare to simulation
p.circle(x = LM_19['Concentration (mM)'][1:]*4, y = LM_19['GFP/OD700'][1:]/150-600)
bokeh.io.show(p)
Changing Parameters to get ttr impact¶
Create Parameter Tabling Functions¶
In [19]:
rx_list = [
'Polymerase Bind to P1d',
'Ribosome binding to ttrS mRNA',
'Ribosome binding to ttrR mRNA',
'ttrS phosphorylation',
'ttR binds to unphosphorylated ttrS',
'ttR binds to phosphorylated ttrS',
'phosphorylation of ttrR',
'dimerization of ttrR~P',
'Dimer bindgint to pTTR',
'Polymerase binding to active pTTR',
'Ribosome binding to GFP mRNA',
]
extra_rx_list = ['transcription', 'translation', 'phosphorylation', 'degradation']
In [20]:
def print_rxn_df(parameters, x0):
params_list_u = [p[0] for p in parameters if p[0].find('u') != -1]
params_list_b = [p[0] for p in parameters if p[0].find('b') != -1]
values_u = [p[1] for p in parameters if p[0].find('u') != -1]
values_b = [p[1] for p in parameters if p[0].find('b') != -1]
params_list_extra = [p[0] for p in parameters if p[0].find('b') == -1 and p[0].find('u') == -1]
values_extra = [p[1] for p in parameters if p[0].find('b') == -1 and p[0].find('u') == -1]
key_list = []
val_list = []
for key, value in x0.items() :
key_list.append(key)
val_list.append(value)
params = pd.DataFrame()
params['reaction'] = rx_list
params['ku'] = values_u
params['kb'] = values_b
extra_params = pd.DataFrame()
extra_params[''] = ['|']*len(values_extra)
extra_params['extra reaction'] = extra_rx_list
extra_params['k'] = values_extra
inits = pd.DataFrame()
inits[''] = ['|']*len(key_list)
inits['initial'] = key_list
inits['val'] = val_list
return pd.concat([params, extra_params, inits], axis=1).fillna("")
Plotting Function¶
In [21]:
def plot_ttr_curves(parameters, x0):
experimental_ttr_concs = np.sort(LM_19['Concentration (mM)'][1:].values)
r=list(bokeh.palettes.viridis(10))
r.reverse()
r = r[3:]
M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
timepoints = np.linspace(0, 1400, 1400)
p = bokeh.plotting.figure(width = 500, height= 400, x_axis_label='Time (hr)', y_axis_label = 'GFP' )
i = 0
for tt_conc in tqdm.tqdm(experimental_ttr_concs):
x0["tt"] = tt_conc #Change my initial condition dictionary
M.set_species(x0)
Results = py_simulate_model(timepoints, Model = M, stochastic = False)
p.line(timepoints/60, Results['GFP'], color = r[i], legend_label = f'{tt_conc}', line_width = 2)
i+=1
p.legend.location = 'top_left'
p.legend.title = 'Tetrathionate mM'
return p
Inverse Tetrathionate Response¶
This set of parameter values gives higher activation when low tetrathionate as opposed to high tetrathionate.
In [22]:
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":250,
"R":20,
}
In [23]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 0.004), ("k2u", 0.1),
("k3b", 0.003), ("k3u", 5),
# Slow reverse reaction on ttrS phosphorylation
("k4b", 0.0005), ("k4u", 0.5),
# ttrR should quickly bind to unphos
("k5b", 5), ("k5u",6),
("k6b", 5), ("k6u", 10),
("k7b", 1), ("k7u", 100),
("k8b", 80), ("k8u", 5),
("k9b", 400), ("k9u", 3),
("k10b", 50), ("k10u",1),
("k11b", 10), ("k11u", 1),
("ktx", .1), ("ktl", 1), ("k_dephos", 0.07), ("delta", .7)]
In [24]:
print_rxn_df(parameters, x0)
Out[24]:
In [25]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
Curves start together, then spread apart¶
In experiments, I found the different induction curves start close together and then fan out towards the end. I attempted find parameters that matched this behavior.
In [26]:
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":300,
"R":10,
}
In [27]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 0.01), ("k2u", 100),
("k3b", 0.3), ("k3u", 500),
# Slow reverse reaction on ttrS phosphorylation
("k4b", 5), ("k4u", 0.5),
# ttrR should quickly bind to unphos
("k5b", 5), ("k5u",6),
("k6b", 0.05), ("k6u", 1),
("k7b", 0.0001), ("k7u", 1),
("k8b", 8), ("k8u", 5),
("k9b", 4), ("k9u", 3),
("k10b", 50), ("k10u",1),
("k11b", 10), ("k11u", 10),
("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .7)]
In [28]:
print_rxn_df(parameters, x0)
Out[28]:
In [29]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
Curves start apart, then come together¶
Now, here is the case where induction curves start together and spread apart.
In [30]:
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":300,
"R":25,
}
In [31]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 0.3), ("k2u", 10),
("k3b", 0.03), ("k3u", 50),
# Slow reverse reaction on ttrS phosphorylation
("k4b", 5), ("k4u", 0.5),
# ttrR should quickly bind to unphos
("k5b", 5), ("k5u",6),
("k6b", 0.05), ("k6u", 1),
("k7b", 0.0001), ("k7u", 1),
("k8b", 8), ("k8u", 5),
("k9b", 4), ("k9u", 3),
("k10b", 500), ("k10u",1),
("k11b", 100), ("k11u", 100),
("ktx", .1), ("ktl", 1), ("k_dephos", 70), ("delta", .7)]
In [32]:
print_rxn_df(parameters, x0)
Out[32]:
In [33]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
Amplified inverse response curve¶
Compared to the previous linear inverse response, the low induction curves grow in a concave curve in this example.
In [34]:
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":300,
"R":1000,
}
In [35]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 0.3), ("k2u", 10),
("k3b", 0.003), ("k3u", 50),
# Slow reverse reaction on ttrS phosphorylation
("k4b", 5), ("k4u", 0.5),
# ttrR should quickly bind to unphos
("k5b", 500), ("k5u",6),
("k6b", 0.05), ("k6u", 1),
("k7b", 0.1), ("k7u", 1),
("k8b", 8), ("k8u", 5),
("k9b", 4), ("k9u", 3),
("k10b", 50), ("k10u",1),
("k11b", 10), ("k11u", 10),
("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .7)]
In [36]:
print_rxn_df(parameters, x0)
Out[36]:
In [37]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
Close match to experiments¶
Now, we choose parameters that give decent quantitative agreement with the experimental data
In [38]:
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":250,
"R":20,
}
In [39]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 3), ("k2u", 10),
("k3b", 0.00003), ("k3u", 5),
# Slow reverse reaction on ttrS phosphorylation
("k4b", 5), ("k4u", 0.5),
# ttrR should quickly bind to unphos
("k5b", 5), ("k5u",6),
("k6b", 5), ("k6u", 1),
("k7b", 0.01), ("k7u", 1),
("k8b", 80), ("k8u", 5),
("k9b", 0.004), ("k9u", 3),
("k10b", 50), ("k10u",1),
("k11b", 10), ("k11u", 1),
("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .9)]
In [40]:
print_rxn_df(parameters, x0)
Out[40]:
In [41]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
RBS Tuning¶
We take the above well-matching plot and perform RBS tuning, like we did experimentally with the ARL pool. We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively.
First, we will define a new plotting function that can take in different parameter values and plot on the same axis
In [42]:
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u):
experimental_ttr_concs = np.sort(LM_19['Concentration (mM)'][1:].values)
r=list(bokeh.palettes.viridis(10))
r.reverse()
r = r[3:]
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1),
("k2b", 3), ("k2u", k2u),
("k3b", 0.00003), ("k3u", k3u),
# Slow reverse reaction on ttrS phosphorylation
("k4b", 5), ("k4u", 0.5),
# ttrR should quickly bind to unphos
("k5b", 5), ("k5u",6),
("k6b", 5), ("k6u", 1),
("k7b", 0.01), ("k7u", 1),
("k8b", 80), ("k8u", 5),
("k9b", 0.004), ("k9u", 3),
("k10b", 50), ("k10u",1),
("k11b", 10), ("k11u", 1),
("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .9)]
M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
timepoints = np.linspace(0, 1400, 1400)
p = bokeh.plotting.figure(width = 500, height= 400, x_axis_label='Time (hr)', y_axis_label = 'GFP', title = f'ttrS (k2u): {k2u}, ttrR (k3u): {k3u}' )
i = 0
for tt_conc in experimental_ttr_concs:
x0["tt"] = tt_conc #Change my initial condition dictionary
M.set_species(x0)
Results = py_simulate_model(timepoints, Model = M, stochastic = False)
p.line(timepoints/60, Results['GFP'], color = r[i], legend_label = f'{tt_conc}', line_width = 2)
i+=1
p.legend.location = 'top_left'
p.legend.title = 'Tetrathionate mM'
p.y_range = Range1d(0,1.9*10**(-8))
return p
In [43]:
# Define higher off rates
k2u = [20, 12, 7]
k3u = [15, 8, 5]
plot_list = []
for i in range(3):
plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i]))
bokeh.io.show(gridplot([plot_list]))
Great, this looks similar to the RBS tuning experiments we did.
Saving and Loading in as SBML¶
I did this as a check that SBML io ran properly.
In [44]:
# #Write the model M to a file
# #Note that bioscrape saves the current state of the model and will not auto-save any changes to it
# file_name = 'SignalTxTl.xml'
# M.write_bioscrape_xml(file_name)
# M_loaded = Model(file_name)
# #M_SBML = import_sbml('sbml_file.xml') #This code won't work without first creating an SBML File
# #And the simulation will look the same as before
# Results = py_simulate_model(timepoints, Model = M_loaded, stochastic = False)
# plt.plot(timepoints, Results["S"]+Results["F1"], label = "S det", color = 'blue')
# plt.plot(timepoints, Results["T"]+Results["T:R"], label = "T det", color = 'cyan')
# plt.plot(timepoints, Results["X"], label = "X det", color = 'purple')
In [45]:
# M.write_sbml_model("TTR_only.xml")
In [46]:
M_ttr = import_sbml('TTR_only.xml')
#Simulate Deterministically and Stochastically
timepoints = np.linspace(0,700,10000)
result_det = py_simulate_model(timepoints, Model = M_ttr)
result_stoch = py_simulate_model(timepoints, Model = M_ttr, stochastic = True)
In [47]:
import matplotlib as mpl
color_list = ['r', 'k', 'b','g','y','m','c']
mpl.rc('axes', prop_cycle=(mpl.cycler('color', color_list) ))
mpl.rc('xtick', labelsize=16)
mpl.rc('ytick', labelsize=16)
plt.figure(figsize = (10, 4))
Out[47]:
In [48]:
#Plot Results
for i in range(len(M_ttr.get_species_list())):
if i < len(color_list):
s = M_ttr.get_species_list()[i]
plt.plot(timepoints, result_det[s], color = color_list[i], label = "Deterministic "+s)
plt.plot(timepoints, result_stoch[s], ":", color = color_list[i], label = "Stochastic "+s)
plt.title('TTR Model')
plt.xlabel('Time', FontSize = 16)
plt.ylabel('Amount', FontSize = 16)
plt.legend()
plt.show()
HRP AND Circuit¶
This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.
In [49]:
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = 'BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("mydna",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_mydna":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
hv.Points(
data=results_melted.loc[results_melted['species'].str.startswith('protein')],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
['species']
).opts(
tools=['hover'],
legend_position = 'right',
width = 800,
#shared_axes=False
).overlay()
Out[49]:
Vary hrpR and hrpS proteins¶
Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.
Get plotting in order¶
Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been
In [50]:
def simulate_heatmap(inducer_1_name,
inducer_1_values,
inducer_2_name,
inducer_2_values,
title,
x0=x0,
normalize=True,
max_to_1 = False,
CRN=CRN_extract_1,
timepoints=timepoints):
GFP_max = pd.DataFrame()
for conc_1 in inducer_1_values:
x0[inducer_1_name] = conc_1
for conc_2 in inducer_2_values:
x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
GFP_max = GFP_max.append({inducer_1_name:conc_1,
inducer_2_name:conc_2,
'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)
data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
columns = inducer_2_name,
values = 'GFP_max')
fig, ax = plt.subplots()
if normalize:
data /= np.max(GFP_max['GFP_max'])
if max_to_1:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.1, vmin = 0, vmax = 1)
else:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.2)
heat_map.set_title(title)
fig.tight_layout()
return fig
Normal Functioning And gate¶
This non-leaky case shows that when we increase protein hrpR and hrpS, we see a continued increase in the amount of GFP produced.
In [51]:
hrpS_values = [0, 0.001, 0.01, 0.1, 0.5, 1, 10]
hrpR_values = [0, 0.0001, 0.001, 0.01, 0.1, 0.5, 1]
fig = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='hrpR, hrpS protein AND gate')
#fig.savefig("test.png")
In [52]:
pd.read_csv('parameters_and.txt', sep="\t", header=0)[49:]
Out[52]:
And Gate with leak in hrpR¶
Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.
In [53]:
pd.read_csv('parameters_and_leak.txt', sep="\t", header=0)[49:]
Out[53]:
In [54]:
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = 'BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("mydna",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "parameters_and_leak.txt")
#print("BioCRNpyler Representation:\n", repr(extract_1_TXTL))
CRN_extract_1 = extract_1_TXTL.compile_crn()
#CRN_model = extract_1_TXTL.get_model()
#print("\nCRN Representation:\n", repr(CRN_extract_1))
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_mydna":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
hv.Points(
data=results_melted.loc[results_melted['species'].str.startswith('protein')],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
['species']
).opts(
tools=['hover'],
legend_position = 'right',
width = 800,
#shared_axes=False
).overlay()
Out[54]:
In [55]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
fig = simulate_heatmap('protein_hrpR',hrpR_values,'protein_hrpS',hrpS_values,x0=x0,normalize=False,title='hrpR, hrpS protein AND gate', CRN=CRN_extract_1)
fig.savefig("leak.png")
And Gate with lots leak in hrpR becomes NAND¶
In [56]:
pd.read_csv('parameters_and_hrpR.txt', sep="\t", header=0)[49:]
Out[56]:
In [57]:
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = 'BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("mydna",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "parameters_and_hrpR.txt")
#print("BioCRNpyler Representation:\n", repr(extract_1_TXTL))
CRN_extract_1 = extract_1_TXTL.compile_crn()
#CRN_model = extract_1_TXTL.get_model()
#print("\nCRN Representation:\n", repr(CRN_extract_1))
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_mydna":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
hv.Points(
data=results_melted.loc[results_melted['species'].str.startswith('protein')],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
['species']
).opts(
tools=['hover'],
legend_position = 'right',
width = 800,
#shared_axes=False
).overlay()
Out[57]:
In [58]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1, 10.0]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1, 10.0]
fig = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=False,title='hrpR, hrpS protein AND gate', CRN=CRN_extract_1)
fig.savefig("nand.png")
AND gate¶
LacI promoter class¶
In [59]:
class LacIPromoter2(Promoter):
def __init__(self, name = "pLac", IPTG = "IPTG",
lacI = "lacI", leak = True, assembly = None,
transcript = None, length = 0, mechanisms = {},
parameters = {}, **keywords):
Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
#RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)
# Transcription is already included in promoter class
self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}
self.leak = leak
self.IPTG = self.set_species(IPTG, material_type = "inducer")
self.lacI = self.set_species(lacI, material_type = "protein")
self.plac_2_lacI = None
self.IPTG_2_lacI = None
def update_species(self):
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
species = []
self.complexes = []
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + pLac <--> 2xlacI:pLac
"""
species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
species += species_b
#print('plac stuff \n',species_b)
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
self.plac_2_lacI = s
if self.leak is not False:
species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
# IPTG can bind to lacI to sequester it
# give citation for 1 being enough
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + IPTG <--> 2xlacI:IPTG
"""
species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
#print('iptg stuff \n', species_b)
species += species_b
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
self.IPTG_2_lacI = s
# Unbound transcription
species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return species
def update_reactions(self):
print('r')
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
self.update_species()
if self.leak != False:
# Will need transcription, plac_2x_lacI, ktx, ku, kb
reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
transcript = self.transcript, protein = self.assembly.protein)
reactions = []
# Repression
# Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
part_id = self.plac_2_lacI.name, cooperativity = 2)
# LacI binding to IPTG
# Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
# .name gets rid of the complex
# .type will give you complex, species, etc
# .attributes gives you properties like phosphorylation state
reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)
'''
Any bit of IPTG destroys the bound structure
IPTG + plac_2xlacI --> IPTG:2xlacI + plac
'''
if isinstance(mech_b, Two_Step_Cooperative_Binding):
kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
else:
kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))
# Unbound transcription plac, kb, ku, ktx
reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return reactions
Sort out param printer¶
In [60]:
def param_printer2(x0, parameters):
inits = pd.DataFrame(list(x0.items()),columns = ['param','value'])
df = pd.DataFrame(list(parameters.items()),columns = ['param','value'])
col_part_id = []
mech = []
for val in df['param']:
if isinstance(val, tuple):
col_part_id.append(val[0])
mech.append(val[1])
else:
col_part_id.append(' ')
mech.append(val)
df['param'] = mech
df['part'] = col_part_id
df.reindex(['part','param','value'],axis=1)
df[''] = ['|']*len(col_part_id)
return pd.concat([df, inits], axis=1).fillna("")
AND simulation¶
In [61]:
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_pttr_hrpR":5,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":10,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":5,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":10,
"protein_Ribo":90.,
"ligand_tt": 100,
"inducer_IPTG": 10}
In [62]:
parameters = {
"kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":2, 'cooperativity':2, #some default parameters
('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0,
('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025, #this is the transcription rate
('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('translation_mm', 'B0030', 'ku'): 10.0, #Unbinding
('translation_mm', 'B0030', 'ktl'): 1.5, #Translation Rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm','kb'):100, #default binding rate!
#('ttrS-P',"kdephos"):1000, #auto dephosphorylation
('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
#('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
#('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
('ttrR-P',"kdephos"):1000, #auto dephosphorylation
('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
('phrpL_hrpS', 'cooperativity'): 1.0,
('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025, #this is the transcription rate
('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('translation_mm', 'B0030', 'ku'): 10.0, #Unbinding
('translation_mm', 'B0030', 'ktl'): 1.5, #Translation Rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm', 'ku'): 100,
('phrpL_hrpR','ku'):100,
('phrpL_hrpS','ku'):100,
('phrpL_hrpS_RNAP', 'ku'): 100.,#we'll say that leak has 10x the unbinding rate
('phrpL_hrpR_RNAP', 'ku'): 100,
('phrpL','ku'):100,
('binding', 'pLac_lacI','ktx'):0.001,
('plac_lacI','kb'):10,
('plac_lacI','ku'):0.001,
}
#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")
# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")
lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")
# # constit_P
# plac_hrpS = DNAassembly(name = "p_hrpS",promoter=RegulatedPromoter("P",["lacI"]),rbs="rbsH",protein="hrpS")
# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")
hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)
CRN_extract_1 = extract_1_TXTL.compile_crn()
In [63]:
param_printer2(x0, parameters)
Out[63]:
In [64]:
#print("\nCRN Representation:\n", CRN_extract_1.pretty_print())
In [65]:
IPTG_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
tt_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
fig = simulate_heatmap('ligand_tt',tt_conc_list,'inducer_IPTG',IPTG_conc_list,x0=x0,normalize=False,title='Compiled AND gate', CRN=CRN_extract_1)
fig.savefig("and_first.png")
Weird nand behavior again...¶
In [66]:
x0 = {"dna_pttr_hrpR":5,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":10,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":5,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":100,
"protein_Ribo":90.,
"ligand_tt": 100,
"inducer_IPTG": 10}
In [67]:
IPTG_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
tt_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
fig = simulate_heatmap('ligand_tt',tt_conc_list,'inducer_IPTG',IPTG_conc_list,x0=x0,normalize=False,title='Compiled AND gate', CRN=CRN_extract_1)
#fig.savefig("and_first.png")
Increasing RNAP, Ribo 10x from initial and¶
In [68]:
x0 = {"dna_pttr_hrpR":5,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":10,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":5,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":100,
"protein_Ribo":900.,
"ligand_tt": 100,
"inducer_IPTG": 10}
In [69]:
param_printer2(x0, parameters)
Out[69]:
In [70]:
IPTG_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
tt_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
fig = simulate_heatmap('ligand_tt',tt_conc_list,'inducer_IPTG',IPTG_conc_list,x0=x0,normalize=False,title='Compiled AND gate', CRN=CRN_extract_1)
#fig.savefig("resource.png")
Computing Environment¶
In [71]:
%load_ext watermark
%watermark -v -p numpy,pandas,scipy,bioscrape,biocrnpyler,bokeh,jupyterlab,holoviews